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Abstract 

We conduct a theoretical study in which we determine the zero-point vacancy concentration in 
solid 4 He at T = OK. To this end, we employ the quantum-classical isomorphism, by which the 
quantum-mechanical probability density function of a system composed of bosons at T = OK can 
be interpreted in terms of a Boltzmann factor of a classical system at finite temperature. By using 
this classical isomorph we apply the methods of classical statistical mechanics to compute the 
vacancy formation free energy and the vacancy concentration in the associated quantum system 
at T = 0. In this context, we focus specifically on the role of anharmonic effects that are expected 
to be non-negligible due to the significant zero-point motion. For this purpose, we compute the 
formation free energies using both the harmonic approximation (HA) as well as reversible-work 
(RW) method, in which all anharmonic effects are taken into account. The results suggest that 
anharmonic effects indeed play a significant role, lowering the classical formation free energy by 
~ 25% and increasing the zero-point vacancy concentration by more than an order of magnitude 
compared to the HA. 

PACS numbers: 61.72.jd, 63.20.Ry, 67.80.B- 
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I. INTRODUCTION 



The unusual properties of the condensed phases of Helium have turned this substance 
into one of the most studied in the literature. An example is the fact that Helium does 
not exhibit a solid phase at ambient pressure, not even at temperatures tending to absolute 
zero. This is due to the weak interactions between the Helium atoms and the large zero-point 
motion. Accordingly, only the application of a significant external pressure, of the order of 
25 atm, allows the formation of a stable solid phase. The high-pressure solid phase of 4 He has 
recently become the subject of intensive investigation after the experimental observations of 
nonclassical rotational inertia (NCRI) 4 -^ 1 ^ and unexpected mechanical behavior.- 

Although 4 He's liquid phase and its display of the phenomenon of superfluidity are com- 
paratively well understood, the comprehension of its solid counterpart involves more subtle 
issues that are not so easily grasped (for a recent review see Ref. 1^). This is reflected, for 
instance, in the fact that the origin of the observed NCRI continues to be a subject of de- 
bate. A possible interpretation of this phenomenon is the existence of a supersolid phase, 
which presents the long-range order and spontaneous translational and rotational symmetry 
breaking characteristic of a solid but displaying superfluid-like behavior. The pioneering 
theoretical proposals of a 4 He supersolid phased predict that vacancies are required for 
quantum crystals to exhibit superfluid behavior. More recently, despite of some contro- 
versy,-* 1 ^ there is evidence 1 ^- that vacancies can indeed exist at zero temperature and that 
they may play a role in the observed NCRI. 

In view of the possible role of lattice vacancies in the unusual behavior of solid 4 He, 
it is of interest to estimate their zero-temperature equilibrium concentrations. Given that 
experimental estimates are not yet available, we need to resort to purely theoretical methods. 
Hodgdon and Stillinger— were the first to conduct such an effort, taking advantage of the 
well-known isomorphism between a quantum-mechanical system composed of Bosons at zero 
temperature and a classical system at finite temperature. This isomorphism is a consequence 
of the fact that the wave function describing the ground state of a system of identical 
Bosons is nodeless.— Accordingly, the associated probability density function can be formally 
written in terms of a classical Boltzmann factor exp(— U/ksT), where U is the classical 
potential energy function, T is the absolute temperature, ks is the Boltzmann constant and 
ksT is to be taken equal to 1. By exploiting this isomorphism it is possible to apply the 
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finite-temperature methods of classical statistical mechanics to estimate the properties of 
a quantum-mechanical Boson system at zero temperature. Specifically, the calculation of 
the zero-point equilibrium vacancy concentration in solid 4 He at zero temperature becomes 
equivalent to the determination of the thermal equilibrium vacancy concentration of the 
associated classical solid at a finite temperature. The latter can be computed by determining 
the finite-temperature vacancy-formation free energy. 

Hodgdon and Stillinger— utilized the standard classical harmonic approximation 
(HA) 1 ^ 1 ^ to determine the vacancy formation free energy and corresponding vacancy con- 
centration. In classical solids, however, the HA is known to be problematic in the presence 
of defects,— which, due to the local breaking of lattice symmetries, enhance the role 
of anharmonic effects. In this context, given the significant zero-point motion in quantum 
solids, such anharmonic effects are also expected to be relevant in the determination of 
the zero-temperature equilibrium concentration of lattice vacancies in solid 4 He. In this 
light, the purpose of the present paper is to evaluate the influence of anharmonic effects on 
the zero-point vacancy concentration in solid 4 He. To this end, following the approach of 
Hodgdon and Stillinger,— we apply the aforementioned isomorphism using an approximate 
description of solid 4 He based on a trial function of the Jastrow-form.— ^ However, in ad- 
dition to applying the HA to compute finite-temperature vacancy-formation free energy in 
the associated classical solid, we also utilize an exact method in which all anharmonic effects 
are explicitly included. 

The remainder of the paper has been organized as follows. Section [III discusses the 
details of the isomorphism between a quantum solid composed of Bosons in its ground state 
and a classical crystal. Furthermore, it gives a description of the computational methods 
used to determine the formation free energy and the corresponding thermal equilibrium 
concentration of vacancies in a classical crystal. Section II I II describes the computational 
details of the calculations and in Section [TV] we discuss the obtained results. We conclude 
with a summary in Section [V] 
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II. METHODOLOGY 



A. The isomorphism between classical and quantum Bose systems 

The simplest trial wave function that is able to describe a system of 4 He atoms is a 
pairwise function of the Jastrow form 

1 - 



*(R) = exp 



(1) 



where R = {ri, r2, . . . , r^} represents the N coordinates of the atoms of the system and u(r) 
is the so-called pseudopotential.— ^ In all our calculations we utilize the particular form 



^ 6 



u(r) = - , (2) 



containing the single parameter b.— Although pseudopotentials that give better variational 
energies for the solid phase do exist,— >2i we have chosen this particular representation be- 
cause it has been already used in HA calculations and for its computational simplicity. 
The methodology utilized to evaluate the influence of anharmonic effects, however, is not 
restricted to this particular form and can be applied to any of the more accurate represen- 
tations proposed in the literature. 

The isomorphism that exists between the probability density of finding a quantum system 
of identical Bosons in the ground-state as described by Eq. (CQ) and a classical crystal of 
potential energy U (R) at temperature T is easily realized by writing 

Z7(R) = 5>(t>,-) , (3) 

i<j 

such that 

W = ap (JW), (4, 



k B T , 

with ksT = 1. Here, and in the remainder of the paper, the potential energy U and 
all quantities that have the dimension of energy are measured in arbitrary energy units. 
Eq. (j3J) can be interpreted as the Boltzmann factor for a configuration R of a classical 
system described by the potential energy function Eq. (151) . For the particular form of the 
pseudopotential used in the present paper, Eq. (j2J), the classical system corresponds to a 
collection of purely repulsive soft spheres. 
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B. Vacancy thermodynamics in classical systems 



The determination of the thermal equilibrium concentration of vacancies in a classical 
crystal requires the minimization of the free energy of a crystal containing N atoms with 
respect to the total number of vacancies n.— Let / be the formation free energy of a mono- 
vacancy in the crystal, defined as 

f = F(N-l)-^F(N), (5) 

where F(N — 1) and F(N) are the Helmholtz free energies of, respectively, a crystal con- 
taining N — 1 atoms and a single vacancy, and a defect-free crystal containing N atoms. 
Accordingly, the free energy of a system consisting of N atoms and n vacancies can be 
written as 



F(N + n) = F(N) +nf- k B T\n 



(N + n)\ 



(6) 



n\N\ 

Note that when the number of vacancies is zero in the above expression, F(N + n) reduces to 
the free energy of the defect-free system. The last term in Eq. (J6j) is just TS con f where S^onf 
represents the configurational entropy associated with the number of ways in which the n 
vacancies can be distributed among the lattice sites. In thermodynamic equilibrium the free 
energy is a minimum and the number of vacancies may be determined by minimizing Eq. (jSJ) 
with respect to n. The resulting thermal equilibrium vacancy concentration c is given by 



n 

c^- = exp 



/ 



k B T 



(7) 



C. Free-energy calculations 



1. Reversible work 



In order to determine the vacancy concentration by means of Eq. (j^J), we need to compute 
the vacancy-formation free energy /. Given that the free energy is a thermal quantity that 
cannot be expressed in terms of an ensemble average,— 1 ^ it cannot be computed directly 
using any Monte Carlo or Molecular Dynamics sampling method. As a result, free energies 
are usually determined using indirect strategies, in which free-energy differences between 
two systems are computed by evaluating the work associated with a reversible process that 
connects these two systems.— This approach has shown to be very effective for the compu- 
tation of defect thermal equilibrium concentrations in classical solids.— Accordingly, by 



virtue of the discussed quantum-classical isomorphism, it should also prove useful for the 
evaluation of zero-point defect concentrations in Boson quantum solids. 

The main idea of the reversible work (RW) method is to take a reference system, for which 
the free energy is known, and connect it to the system of interest by means of a coupling 
parameter A. A typical functional form of this coupling is given by the Hamiltonian^i 

H(X) = XH + (1 - A)# ref , (8) 

where Hq and H rc f represent the Hamiltonians of the system of interest and of reference, 
respectively. Note that this form allows a continuous switching between if re f and H by 
varying the parameter A, varied between and 1. The free-energy difference between the 
system of interest and the reference is then given by the reversible work W re ^ 

W rev = F - F TC{ = J dX /|^\ , (9) 

where the brackets indicate an equilibrium average in the relevant statistical ensemble (i.e., 
canonical, isobaric-isothermal, etc.). The integration represents the calculation of the to- 
tal work done by the generalized force dH/dX. Since the integration involves equilibrium 
averages of the system at all times, it reflects a reversible process. 

In principle, the numerical evaluation of the work integral Eq. (Q requires a series of 
independent equilibrium simulations, each carried out at a different value of the coupling 
parameter A between and 1.— In practice, however, it has shown beneficial to estimate 
the work integral along a single, nonequilibrium simulation during which the value of A 
changes dynamically. In this fashion, the reversible work W rev is estimated in terms of the 
irreversible work estimator 

'dH({ri},xy 



r ^sim 


~dX~ 




/ d£ 




Jo 


dt_ 


t 1 . 



OX 



(10) 

J A(t') 



where t' represents the time coordinate that describes the dynamical evolution of the coupling 
parameter X(t), and t S i m is the total duration of the switching process. 

Given that the dynamical process above is intrinsically irreversible, the irreversible work 
estimator (fTUj) is biased, subject to a positive systematic error that is associated with the 
dissipative entropy production inherent to nonequilibrium processes.— i22i2I As a result, it 
will represent an upper bound to the value of the reversible work W Tev and, consequently, 
the free-energy difference F — F TC f. Fortunately, however, the systematic error can be 
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readily eliminated, as long as the nonequilibrium process remains within the regime of 
linear response. In this regime, the systematic error is independent of the direction of the 
switching process, i.e., the entropy production is equal for the forward (A = — > 1) and 
backward (A = 1 — > 0), processes.— 1 ^ In this fashion, we can obtain an unbiased estimate 
for W Tev according to 

1 



[W irr (A = - 1) - WWA = 1-0)] 



(11) 



which is subject to statistical errors only.— ^ 

The initial step toward computing the formation free energy of a vacancy, Eq. (jSJ), uti- 
lizing the RW approach involves the calculation of the free energy of a defect-free crystal 
containing N atoms of 4 He, F(N). The strategy will be to use the classical Einstein crystal, 
of which the free energy FEinst is known analytically, as a reference system.— In this manner, 
using the classical potential energy expression of Eq. ([3]), the particular form of the coupled 
Hamiltonian (EI) becomes 



H (R; A) = A 



' N 

£ 

,i<j 



A) 



i N 

Z/(R (0) ) + ±I> (r« 

1=1 



,(0) 



(12) 



where represents the equilibrium lattice positions of N atoms, U(R^) is the static 
contribution to the potential energy and «j is a spring constant. 

Since W Tev is the difference between the Helmholtz free energy of the system of interest 
and the reference system, the free energy of the defect-free crystal is given by 



F(N) = W„ 



Einst 



(13) 



It is convenient to fix the center of mass of the system during the switching process.— In 
the presence of this constraint, the free energy of the Einstein crystal is given by 

N 



Fm 



Einst 



W(R 



(0)> 



°-k B T^ 



8=1 



2irk B T 



knT In 

2 



_8vr 2 (fc B T) 2 _ ' 



(14) 



The value of the spring constants Ki is adjusted such that the mean-square displacement 
of the harmonic oscillator is approximately equal that of the particles in the system of 
interest .— 

In the second step we estimate the free energy F(N — 1) of the system containing a single 
vacancy. In order to apply the RW method for this purpose, we use a strategy that is similar 



to the one employed to compute F(N). In this case, the parameter A controls the strength 
of the interaction between a single particle and all the others. The coupled Hamiltonian for 
this process is given by 



where the index k stands for the particle that is turned on and off as A is varied between 
and 1. The final term describes a single harmonic oscillator, which is the final state 
of particle k when it has been totally decoupled from the remainder of the system. The 
coupling to the harmonic oscillator is implemented to prevent the decoupled particle k from 
drifting freely through the system.— As before, the center of mass is held fixed during the 
switching simulations. Furthermore, it is convenient^ to maintain the particle density p of 
the system constant. This is done by varying the size L of the simulation box as the value 
of A is changed. Specifically, we set 



such that the density of interacting particles is equal to p both for A = and A = 1. 

A final concern involves the applicability of the approach at elevated temperatures, where 
vacancy diffusion may interfere with the RW process. In this situation, one of the twelve 
nearest-neighbor atoms of the decoupled particle may move into the vacant lattice site by 
means of a diffusion event, which leads to singularities in the driving force.— An adequate 
solution to this problem is to constrain the motion of the neighboring atoms such that they 
are forced to remain within their respective Wigner-Seitz primitive cells.— 

From the coupled Hamiltonian in Eq. (jl5p we determine the driving force dH/dX and 
compute the reversible work W r ' ev required to turn off the interactions between particle k 
and the remainder of the system and turn it into an independent harmonic oscillator. As 
in the previous case, the switching process is carried out in both directions to eliminate the 
systematic error. Accordingly, F(N — 1) is obtained through 



N , , v 6 N / , x 6 1 x9 

*(M;AH£(A) + a E (±) +(1 _ A) i B ( r4 _ r2 ») , (15) 



(i^k) (i^k) 




(16) 



F(N - 1) = W; cv + F(N) - F c 



OSC ) 



(17) 



where F osc is the free energy of a single harmonic oscillator, given by— 




(18) 
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Finally, combining the results of Eqs. (fl3|) . and (ITTj) we compute the vacancy formation 
free energy through Eq. (jSJ). The corresponding zero-point vacancy concentration is then 
given by Eq. ([7j). 

2. Harmonic approximation 

The harmonic approximation^ is a method that allows one to estimate the vacancy- 
formation free energy without sampling configurations from a statistical ensemble. Hodgdon 
and Stillinger— used this approach to estimate the concentration of point defects in solid 
4 He described by a Jastrow wave function. The central idea of this approximation is to write 
the free energy as the sum of static potential energy plus a finite-temperature part involving 
the vibrational normal mode frequencies o>j, where i — 1, . . . , 3N — 3. 

In this manner,— the vacancy formation / is given by 

37V-6 

+ ksT In 

i=l 

where U(N — 1) represent the static potential energies of the crystal containing a single 
vacancy. Similarly, the last two terms represent the vibrational parts of the crystals with 
and without vacancy, respectively. In practice the vibrational frequencies are obtained by 
diagonalizing the Hessian matrix of the potential-energy functions.— 

III. COMPUTATIONAL DETAILS 

All calculations were carried out using fee cubic computational cells containing numbers 
of particles N varying between 108 and 1372. To eliminate surface effects, we apply standard 
periodic boundary conditions and invoke the minimum-image convention.— In order to be 
consistent with a continuous and differentiable wave function^ we force the pseudopotential 
to be equal to zero for interatomic distances r = L/2. A convenient way to do is to slightly 
modifying the pseudopotential u(r) of Eq. (j2J) according to 

u(r) -> u(r) = u{r) + u{L - r) - 2u(L/2) . (20) 



Ui{N 



N -2 
N- 1 



UbT In 



k B T 



(19) 
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In most of our calculations we have chosen b = 1.72a, with a = 2.556A, for which the classical 
system crystallizes into a fee structure, as determined by Prestipino et al..— In practice, 
the system remains a metastable fee solid for smaller values of b, and we have done a few 
additional calculations for smaller 6-values. Hodgdon and Stillinger,— for instance, chose 
the parameter b = 1.67<j for which the OK melting density of 4 He equals the melting density 
of the classical crystal as obtained by Hoover et al..— We find, however, that this value of b 
does not describe a metastable fee solid and leads to a liquid phase.— 

We apply the Metropolis algorithm to sample the configurations from |\l/(i?)| 2 . In the 
RW switching simulations, the switching parameter is varied linearly according to A = t/t S i m 
(forward process) or A = 1 —t/t S i m (backward process). Here, t is an index of the MC sweep 
in the process and t S j m represents the total number of sweeps of the switching process. It 
is important to remember that one must always start an RW switching process from an 
equilibrated system. Accordingly, every RW switching run, for both switching directions, 
was preceded by an equilibration run of at least 5 x 10 3 MC sweeps before computing the 
work Eq. (fTUl) . In order to guarantee that the RW results obtained from nonequilibrium pro- 
cesses are bias-free, we compute the irreversible work estimators Eq. ( flOl) for both switching 
directions as a function of the total switching duration t B i m and monitor the convergence of 
Eq. pip . For each value of t s[m we determine averages of the irreversible work estimator 
over 30 switching processes. 

IV. RESULTS AND DISCUSSION 

Fig. [1] shows the results of the RW calculations for the defect-free crystal (panel a) as well 
as for the reversible vacancy insertion process (panel b) obtained for the cell containing 500 
atoms, with ksT = 1- The graphs display the average values of the forward and backward 
irreversible work estimators, as well as the unbiased estimator Eq. ( ITTj) . as a function of the 
total process duration t sim . The results show that t S i m ~ 10 4 sweeps is sufficiently long for 
the linear response regime to be reached in both cases, for which the unbiased estimator 
reaches convergence. The horizontal lines on the graphs are averages of the last four values 
of the reversible work in that region. Based on these results, all subsequent RW switching 
simulations based on different numbers of particles and temperatures were carried out with 
t sim = 1.5 x 10 4 MC sweeps. 
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The results depicted in Fig. [2] show a comparison between the vacancy-formation free 
energy as a function of temperature as obtained using the RW and HA methods, for a cell 
containing 500 particles. The data clearly demonstrate the significance of anharmonicity. 
As expected, for low temperatures, anharmonic effects are small and both the HA and 
the RW method converge to the same formation free-energy value for zero temperature. 
As temperature increases, however, anharmonicity becomes increasingly important and the 
deviation between the RW and HA results becomes significant. Specifically, for ksT = 1, 
for which the isomorphism between the classical and quantum description of the system 
is valid, anharmonic effects are responsible for a reduction of about 25% in the value of 
the vacancy- formation free energy compared to the HA result. The anharmonic effects are 
associated with the significant zero-point motion in the quantum system at T = K. The 
extent of this motion can be approximately estimate in terms of the formation entropy of 
the classical system, s = — at ksT = 1, which can be computed from the slope of the 
formation free-energy curves. The resulting values give s = l.lks and s = 6.0/cb for the 
HA and RW results, respectively. The difference of almost a factor of six is another clear 
indication of the importance of anharmonicity in the system. 

To assess the influence of finite-size effects, we computed the RW vacancy-formation free 
energy at ksT = 1 for different cell sizes, ranging from N = 108 to N = 1372 atoms. The 
results are shown in Fig. [31 that plots the formation free energy as a function of 1/N. The 
inset shows the full formation free-energy versus temperature curves for different cell sizes. 
The results allow us to extrapolate the vacancy-formation free energy to 1/N — > 0, which 
gives / = 12.24 for infinite cell size. 

The determination of the extrapolated vacancy-formation free energy at k^T = 1 now 
allows us to compute the zero-point vacancy concentration for the quantum system de- 
scribed by Eqs. ([I]) and (|2J) by using Eq. Aside from the RW data, we have also used 
extrapolated HA formation free-energy data (not shown in Fig. [3]) for the calculation of 
the equilibrium concentration. In addition, we have also considered the dependence of the 
vacancy concentration on the value of the parameter b. The results are shown in Fig. HI 
which plots the extrapolated vacancy-formation free energy as well as the associated vacancy 
concentration as a function of the parameter b, using both the RW and HA methods. The 
effects of anharmonicity are clearly visible, reducing the formation free energy by ~ 20 — 25% 
and increasing the vacancy concentration by more than an order of magnitude. The results 
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are also seen to be quite sensitive to the value of the parameter b, which is to be expected 
given that variations in b can be interpreted in terms of variations in the temperature of the 
associated classical isomorph. 

We have also determined the vacancy concentration for densities above the A He melting 
density. For this purpose we have computed the extrapolated vacancy-formation free energies 
as a function of the density of the classical soft-sphere isomorph at ksT = 1 and for b = 
1.72a. The results are shown in Figure [5j It is clear that c decreases strongly with increasing 
density. Specifically, a ~ 20% density increase with respect to the melting density p me it = 
0.479a" 3 leads to a vacancy concentration reduction of 4-5 orders of magnitude. As expected, 
an increasing density also leads to a slight reduction of the anharmonic effects. 

As a final point, we have monitored the frequency of diffusion attempts of the twelve 
nearest-neighbors in the RW vacancy creation simulations. To this end, for b = 1.72<r, we 
recorded the number of times a trial move of any of the 12 nearest neighbors was rejected 
because of a Wigner-Seitz (WS) cell boundary crossing event as a function of temperature 
and the number of particles N in the cell. The results are summarized in Fig. El which 
presents an Arrhenius plot of the rejection probability Pws as a function of temperature. 
All the data fits very well to a straight line, irrespective of the number of atoms in the cell, 
allowing us to estimate the energy barrier AE that needs to be overcome for a vacancy to 
diffuse to a neighboring lattice site in the classical isomorph. Using Pws ~ exp(— AP/Zc^T), 
the slope of the Arrhenius plot gives AE = 2.21. While knowledge of the height of the 
energy barrier permits an estimation of the jump frequency in the classical isomorph, this 
evidently is not the case for the quantum system, for which the determination of tunneling 
frequencies requires knowledge of the shape of the potential-energy barrier. 

V. CONCLUSIONS 

In this paper we demonstrate the importance of anharmonic effects in calculating the 
zero-point vacancy concentration in solid 4 He using the RW and HA methods based on the 
well-known quantum/classical isomorphism for the ground state of a many-body system of 
bosons. We find that anharmonicity leads to a decrease of the formation free energy of about 
25%, resulting in a vacancy concentration that is more than an order of magnitude larger 
compared to the harmonic results. 
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We conclude by noting that the RW methodology for the calculation of free energies in 
classical isomorph is not restricted to the simple pseudopotential Jastrow form employed in 
this paper. In addition it can also be used to investigate the concentration of other species 
of defects such as interstitials and interactions between them. 
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FIG. 1: Irreversible work calculations associated with the process of transforming the interacting, 
defect-free solid into an Einstein crystal (Panel a) and the process of introducing a monovacancy 
into the defect-free solid (Panel b) as a function of the process duration i s j m measured in MC 
sweeps. Results are displayed for the forward processes (A), backward processes (0)> as weu as 
the corresponding unbiased estimator of Eq. (jlip (*). Each point was obtained as the average 
over thirty independent simulations carried out N = 500 particles. The solid line is an average 
of the last four values of the unbiased work estimator, for which the processes have reached the 
linear-response regime. 




FIG. 2: Vacancy-formation free energy as a function of temperature. Results were obtained for 
b = 1.72cr and N = 500. HA (x) and RW (Q) results. The lines are guides to the eyes. 
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FIG. 3: Vacancy-formation free energy as a function of system size at ksT = 1. The straight 
was obtained through a least-square linear fit. Results obtained for b = 1.72a and N = 108 (□), 
N = 256 (■), TV = 500 (O), N = 864 (A) and N = 1372 (y). Inset shows the formation free 
energy as a function of k^T for the mentioned system sizes. 
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FIG. 4: Extrapolated vacancy-formation free energy as a function of b at &bT = 1.0. Results 
obtained using the RW method (A) and the HA (x). Lines are plotted to guide the eyes. 
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FIG. 5: Concentration of vacancies as a function of density for b = 1.72a. Results obtained using 
RW method (a) and HA(x). Lines are plotted to guide the eyes. 
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FIG. 6: Arrhenius plot of Wigner-Seitz rejection probability as a function of temperature. Results 
were obtained for b = 1.72a and N = 108 (□), N = 256 (■), N = 500 (O), N = 864 (A) and 
N = 1372 (v)- The line represents the linear regression to the results obtained with N = 1372. 
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